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Imaginary time Gaussian dynamics of the A13 cluster 
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Semiclassical Gaussian approximations to the Boltzmann operator have become an important tool for the 
investigation of thermodynamic properties of clusters of atoms at low temperatures. Usually, numerically 
expensive thawed Gaussian variants are applied. In this article, we introduce a numerically much cheaper 
frozen Gaussian approximation to the imaginary time propagator with a width matrix especially suited for the 
dynamics of clusters. The quality of the results is comparable to that of thawed Gaussian methods based on 
the single-particle ansatz. We apply the method to the argon trimer and investigate the dissociation process 
of the cluster. The results clearly show a classical-like transition from a bounded moiety to three free particles 
at a temperature T ~ 20 K, whereas previous studies of the system were not able to resolve this transition. 
Quantum effects, i.e., differences with the purely classical case manifest themselves in the low-temperature 
behavior of the mean energy and specific heat as well as in a slight shift of the transition temperature. We 
also discuss the influence of an artificial confinement of the atoms usually introduced to converge numerical 
computations. The results show that restrictive confinements often implemented in studies of clusters can 
influence the thermodynamic properties drastically. This finding may have implications on other studies of 
atomic clusters. 

PACS numbers: 36.40.-c, 03.65. Sq, 05.30.-d 



I. INTRODUCTION 



Rare gas atomic clusters are a topic of ongoing re- 
search partially due to the rich variety of their thermo- 
dynamic properties. Extensive studies have been carried 
out on structural transformations or phase transitions^— . 
Properties of particular interest include the mean en- 
ergy and specific heat. Computations on clusters of light 
atoms, e.g., Nei3 and Ne3gi»2, predict novel low tempera- 
ture quantum effects such as liquid-like zero temperature 
structures of Ne38 as compared to a solid-like structure 
predicted from classical mechanics^. 

The features in the mean energies or specific heats of 
such systems appear usually at low temperatures so that 
accurate quantum mechanical computational methods 
are essential. Accurate calculations for multidimensional 
systems, however, are challenging. Path-integral Monte 
Carlo methods^— have been used to investigate rare gas 
clusters with up to a few dozen atoms, however, they 
become expensive for low temperatures so that approx- 
imations are necessary. Recently, new variants of semi- 
classical initial value representations have been adopted 
to problems involving the Boltzmann (imaginary time) 
operator exp(—/3H). The time evolved Gaussian method 
developed by Mandelshtam and coworkers 3 ^ has suc- 
cessfully been applied to atomic clusters^— i 12 i 13 and dis- 
sipative systems^. It is based on the imaginary time 



propagation of a Gaussian wave packet of the form 

(x\g) = (^ N \detG(r)\)- 1/4 
xexpf-i[ a; -q(r)] T G(T)- 1 [ a; -q(T)]+7(r)) , (1) 
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where for a cluster with iV atoms the vectors are 3N- 
dimensional and G(r) is a 3N x 37V-dimensional sym- 
metric matrix of width parameters. 

The time evolved Gaussian method can also become 
expensive for high dimensional systems. It belongs to 
the so-called thawed Gaussian methods, where the ma- 
trix of Gaussian width parameters G(r) changes with 
time. The number of the resulting equations of motion 
scales with N 2 . This can become difficult when dealing 
with clusters of several dozen atoms. Thus, Mandelsh- 
tam and coworkers introduced the single-particle ansatz 3 , 
in which the width matrix G(r) is reduced to a block- 
diagonal structure by only taking into account correla- 
tions of the coordinates of one particle and ignoring the 
inter-particle connections. With this approximation one 
has linear scaling (6-/V) with the number of atoms. 

From a numerical point of view it is much cheaper to 
use frozen Gaussian representations of the thermal op- 
erator, in which the time-dependent width matrix G(r) 
in Eq. (fTJ) is replaced by a constant matrix. The num- 
ber of the remaining equations of motion which have to 
be solved scale as 37V. Such a formalism has been pro- 
vided by Zhang et al?^. One objective of this article is 
to apply the frozen Gaussian approach to the Boltzmann 
operator for a cluster of atoms. We will show that by an 
adequate nondiagonal choice of the constant width ma- 
trix, the frozen Gaussian method can accurately describe 
the mean energy, the specific heat, and signatures of dis- 
sociation processes. To do so, we will present a simple 
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procedure to find a well suited shape for the width ma- 
trix. The results we obtain are of the same quality as 
the single-particle ansatz thawed Gaussian methodology, 
even though the frozen Gaussian variant leaves much less 
freedom to the Gaussian wave packet (constant width) . 

One cluster which has attracted the interest of theoret- 
ical investigations for a long time is the argon trimer— ~— , 
whose dissociation process has been discussed very re- 
cently in an extensive study 6 . In spite of its apparent 
simplicity with only three atoms involved, the thermody- 
namic properties at low temperatures T < 40 K still in- 
clude open questions. In particular, path-integral Monte 
Carlo calculations of the system^ indicate a dissociation 
of the three atoms at temperatures T <J 35 K but cannot 
distinguish this process from structural changes. With 
the semiclassical Gaussian approximations discussed in 
this article we are able to provide well converged numer- 
ical mean energies and specific heats which exhibit an 
unambiguous classical-like dissociation at T « 20 K. We 
will also address the question of how important quan- 
tum effects are for the dissociation. Influences on the 
low temperature behavior of the mean energy and the 
specific heat will become observable and we will see that 
the transition temperature is shifted to a slightly lower 
value. 

One focus of our discussion will be on the influence of 
an artificial confinement of the atoms. For converging nu- 
merical computations of the thermodynamic properties it 
is often necessary to restrict the configuration space to 
a certain volume by introducing an additional confining 
potential^, which in practical applications usually is cho- 
sen to be very restrictivei^^^. In this article we will 
show that such a restrictive choice can have a drastic 
influence on the dissociation process as was already dis- 
cussed in a classical context many years ago^. This may 
have implications on other studies of atomic clusters. 

The article is organized as follows. In Sec. [TT| we in- 
troduce the Gaussian semiclassical approximation to the 
thermal operator. We review the thawed Gaussian (Sec. 
Ill A[) propagator formalism and develop a new multidi- 
mensional form of its frozen Gaussian counterpart (Sec. 
Ill Bp capable of competing with thawed Gaussian meth- 
ods. The results for the argon trimer based on these semi- 
classical Gaussian methods are then presented in Sec. Mil 
After introducing the system (Sec. IIII A|) and comparing 
the Gaussian methods (Sec. IIII Bp we discuss the influ- 
ence of the confining potential on the thermodynamic 
properties (Sec. IIII CI) and investigate the dissociation in 
the classical and the quantum case ([III Dp . Conclusions 
are drawn in Sec. IIVI 



II. THERMAL OPERATOR FOR CLUSTERS AND 
GAUSSIAN APPROXIMATIONS 

We consider a cluster of N atoms with only internal 
forces depending on the distance between the atoms, i.e., 



the Hamiltonian in mass scaled coordinates has the form 

fc2 N 

^-yE^ + E^-^l)' ( 2 ) 

i— 1 j<i 

where Aj is the Laplacian of particle i and V(|rj— Tj\) de- 
scribes the two-body interaction between particles i and 
j, whose positions are given by the vectors r;. In the 
Gaussian representations used in this article it is neces- 
sary to evaluate integrals over a product of the potential 
with a coherent state, which are typically of the form 

/oo 
dx 3N (x\g({ yi })) 2 h(x), (3) 
-oo 

where (x\g({yi})) is a normalized coherent state in x, 
which depends on a set of parameters {y^} usually in- 
cluding Gaussian positions q and a width matrix G, e.g., 

(x\g ({ Vl } = {q,G})) = (7t 3N \det G\y lfi 

xexp^x-qfG-^x-q^ . (4) 

The function h(x) stands for the potential or one of its 
derivatives. The 3iV-dimensional vectors x and q com- 
bine the coordinates of all N atoms. It is essential for 
practical applications to reduce the numerical integra- 
tions as much as possible. Based on the facts that any 
central potential can be fitted by a sum of Gaussians and 
that a Gaussian in the distance = jr, — Tj\ centered at 
the origin remains a Gaussian in Cartesian coordinates, 
Frantsuzov et al£ suggested the implementation of the 
interaction potential in terms of sums of Gaussians, 

V{\n - r 3 \) = c P e~ a ^ , r iS = \n - r 3 1, (5) 
p 

so that all integrals of the form ^ can be evaluated ana- 
lytically. This is of great advantage in numerical compu- 
tations and is used for the work presented in this article. 

To investigate the thermodynamic properties of the 
cluster we calculate the partition function Z(f3) by eval- 
uating Gaussian initial value representations of the ther- 
mal operator 

K{0) = e^ H , (6) 

where j3 = l/(kT) is the inverse temperature and Z(J5) = 
Tr(K(j3)). We are interested in the mean energy E — 
kT 2 d\nZ/dT and the specific heat C = dE/dT. In our 
calculations we use two different semiclassical propaga- 
tors based on a frozen and on a thawed Gaussian repre- 
sentation, where in both cases the Bloch equation 

- -^\q ,r) = H\q ,T) (7) 

connected with the propagator ([5]) is approximately 
solved for a coherent state |qr , t) sa |p({j/i},r)) with ei- 
ther constant or variable Gaussian width parameters. 
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A. Thawed Gaussian representation 

The thawed Gaussian representation of the thermal op- 
erator is the most versatile since it allows both the po- 
sitions and widths of the Gaussian wave packet to vary 
with time. We consider the symmetrized time evolved 
Gaussian approximation (TEGA) suggested by Frantsu- 
zov et al.— 



{x\K tg {t)\x') 



dq 3N exp[2 7 (r/2)] 
(27r) 3JV det[G(r/2)] 



X exp ( --[* - q(r/2)} T G(r/2)- 1 [x - q(r/2)] 



exp ( --[*' - q(r/2)] T G(r/2)- 1 [ : r' - q(r/2)] ) , 



(8) 



which is constructed from the coherent state 



(x\g({Vi} = {g(r),G(r)})) = (x\g(q(r), G(r))} 



(7r 3W |detG(r)|) 



-1/4 



x e xp\--[x-q(T)] l G(T)- 1 [x-q(T)] j . (9) 

One can then readily write down the partition function 
as: 



dq 3N exp[2 7 (r/2)] 



(10) 



The width matrix G(r) is symmetric positive definite. 
The Gaussian parameters follow the equations of motion 
in imaginary time r, 



dr 
d 



G(r) = -G(t)(VV t 1/( 9 (t)))G(t) + h 2 l, (11a) 



h q(r) = -G(r)(VV(q(r))), (lib) 
d 7(r) = -^Tr [(VV T V(q(r)))G(r)] (V(q(r))), 



dT 



(lie) 



where (...) represents Gaussian averaged quantities of 
the form ([3]), which can be evaluated analytically for a 
potential ([5]) expressed in terms of Gaussians 3 -, and 1 is 
the 3iV x 3-/V-dimensional identity matrix. The boundary 
conditions 



q{r w 0) = q , 

7 (t « 0) = -y(qf )r, 



G(r ps 0) = fi 2 lr, 



(12) 



are derived by demanding that in the limit r — > the 
Gaussian approximation reduces to the identity operator. 

In the framework of a Gaussian propagator the thawed 
Gaussian representation is usually the most accurate ap- 
proximation to the exact quantum result due to the large 
freedom in the parameters, as has recently been demon- 
strated for a double well potential 2 ^. However, it is also 



the numerically most expensive method. The number 
of equations of motion for the width matrix (|lla[) scales 
with TV 2 , and the matrix operations at each time step 
even scale with iV 3 . This drastic increase in the required 
computing resources is the most critical drawback of the 
method. An attempt for combining the advantages of 
a thawed Gaussian propagator, where some matrix el- 
ements are still governed by the equations of motion 
(|llap - (|llcp . and avoiding the drawback of the expensive 
numerical effort to evaluate it, is achieved with the so- 
called "single-particle ansatz" of Frantsuzov et al* . This 
ansatz, or variations of it, have been applied to several 
types of clustersi^ii 3 -. It uses a block-diagonal matrix 
G(r), where 3x3 symmetric matrices representing one 
particle along the diagonal are the only non-vanishing 
matrix elements. Then the equations of motion (|lla[) - 
(jllcp are only solved for the 3x3 blocks and only the cor- 
responding 3x3 blocks of (VV 1 V(q(T))) are included. In 
the single-particle ansatz the number of equations scales 
with N instead of TV 2 , however, one loses information in 
the non-diagonal 3x3 blocks, which are set to 0. Since 
the Gaussian propagators are in practical applications 
usually evaluated in Cartesian coordinates, in which the 
motions of the particles do not separate, important cor- 
relations between the particles are ignored. Thus, one 
expects that compared to the case of a full matrix the 
single-particle ansatz may lead to results of poorer qual- 
ity. 

In what follows we will call the full matrix variant 
of the thawed Gaussian propagator FC-TG (fully cou- 
pled thawed Gaussian, also referred to as the "fully 
coupled variational-Gaussian-wave-packet Monte Carlo" 
in Ref. l3J) and the single-particle ansatz will be called 
SP-TG (single-particle thawed Gaussian, "single-particle 
variational-Gaussian-wave-packet Monte Carlo" in Ref. 
3). For the argon trimer we will compare these respec- 
tive approximations for the thermodynamic properties 
derived from the partition function with two variants of 
a frozen Gaussian propagator. 



B. Frozen Gaussian representation 

The frozen Gaussian representation of the thermal op- 
erator suggested by Zhang et al.— is based on a multidi- 
mensional frozen Gaussian coherent state 

(x\g({ Vi } = {p(r),q(r),r})> = (x\g(p(r), q(r),T)) 

det(r)\ 1/4 / 1. , , lT „ . 
— ^ij exp(--[x-q(r)} T T[x-q(T)} 



-p^r)-[x-q(r)}), (13) 



where T is in general a 3iV x 3./V-dimensional constant 
width matrix with positive eigenvalues, and q(r) and 
p(r) describe the dynamical variables. The symmetrized 
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frozen Gaussian approximation to the propagator reads 
(x'\K fg (t)\x) = det(r)exp ^-^Tr(r)r^) 



det (2 [1 - exp(-ft 2 rr)]" 



\ oxp ( -i[jc' - x] T r[tanh(ft 2 rr/2)] _x [a;' - x] 

r °° da 3N ( r /2 
I ^L—expl-2 / dT<V(g(r))> 



\3N 



oo 



-[x-q(r/2)fr[x-q(r/2)}) (14) 



with S = (x' + x)/2 and the Gaussian averaged potential 
(V(g(r))) is as defined in Eq. (|3J). Taking the trace yields 
the partition function^ 



Z fg {t) = Tr [^fg(t)] = Vdct(r) cxp ( -— Tr(r)r 



det([l-exp(-ft 2 rr)]" 



r °° do 3Ar / / >t/2 \ 

x/ iexp -2/ dr(y(q(r))) . (15) 



OO 



(2^)^/2 



The numerical evaluation is relatively simple since one 
only needs to solve the 3N imaginary time equations of 
motion 



dq(r) 
dr 



-r- l <W(g(r))) 



(16) 



for the Gaussian positions q(r) and only one configura- 
tion space integration over the initial positions q(j = 0) 
has to be performed. 

As in the case of the SP-TG propagator the numerical 
scaling in the evaluation of the two-body potential terms 
is A 2 . The numerical advantage of the frozen Gaussian 
propagator as compared to the FC-TG or SP-TG meth- 
ods is due to the constant width matrix, i.e., one has only 
to propagate the equations of motion (|16l) . whose num- 
ber scales with N. Additionally, functions of the width 
matrix T can be evaluated in advance and do not have 
to be repeated at every time step since T does not evolve 
in time. On the other hand, the width matrix T is a pa- 
rameter of the system and its actual choice has a critical 
impact on the quality of the results^. It is not trivial to 
find a good choice of T, however, the problem simplifies 
when all particles are identical. 

For N identical particles the simplest structure for T 
is a diagonal matrix with identical width elements, i.e., 
a multiple of the 37V x 3N identity matrix, 



ri = n, 



(17) 



with only one parameter T. This ansatz treats all par- 
ticles equally, is very simple to implement, and has the 



lowest numerical cost due to the diagonal structure of 
the matrix in Cartesian coordinates. However, it ignores 
the fact that a correct description of the cluster has to 
contain both the free motion of the center of mass and 
the relative motion determined by the particle-particle 
interaction potential ([5]). 

In a frozen Gaussian approximation the exact parti- 
tion function of a free particle is obtained in the limit of 
a vanishing Gaussian width [c£, e.g., Eq. (|15p], whereas 
the optimum width for the relative coordinates can be de- 
duced from a harmonic approximation around the min- 
imum of the particle-particle interaction potential and 
has a finite value. A separation of the free motion of 
the center of mass from the internal degrees of freedom 
should be avoided since it complicates the structure of the 
equations of motion by introducing numerically more ex- 
pensive terms such as a non-diagonal mass matrix. Con- 
sidering the thawed Gaussian propagators we note that 
the FC-TG is capable of correctly describing the free cen- 
ter of mass motion when it is combined with an internal 
potential independently of the choice of the coordinate 
system, whereas this is not fulfilled for the SP-TG^ 2 -. 

In the following we suggest a procedure based on an 
adequate choice of the width matrix which allows for a 
correct description of the free center of mass motion with- 
out changing the structure of the equations. To simplify, 
we restrict our description to the case of three particles 
relevant to this article, a generalization to an arbitrary 
number of particles is straightforward. It is plausible that 
in a system of coordinates Ri for the center of mass 



-Rem = ^ (n + r 2 + r 3 ) 



and the two relative positions 



Ri 
R2 



ri 



T2, 

r 3 



(18a) 



(18b) 
(18c) 



a diagonal matrix structure is a good choice. In these 
coordinates the center of mass is separated and we intro- 
duce the Gaussian width parameter D\ for its motion. 
The Gaussian approximation becomes exact for the cen- 
ter of mass motion in the limit D\ — > 0. Since all parti- 
cles are equal and there is no motivation for distinguish- 
ing between the propagation of the individual relative 
coordinates, we use one single parameter D 2 for the re- 
maining coordinates. The matrix which then is applied 
to the coordinates (|18a|) - ()18c|) is 



\Di 
- , D 2 
D 2/ 



(19) 



where D\ and D 2 are 3x3 diagonal matrices with coef- 
ficients D\ and D 2 , respectively, and is a 3 x 3 matrix 
of zeros. The most efficient way to evaluate the frozen 
Gaussian partition function is to keep its structure (|15[) 
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TABLE I. Parameters used in the Gaussian fit ([5]) of the 
Morse potential (f2Tj) . 



V 


c p [cm -1 ] 


a p [A ] 


1 


3.296 x 10 s 


0.6551 


2 


-1.279 x 10 3 


0.1616 


3 


-9.946 x 10 3 


6.0600 



in Cartesian coordinates and to transform the width ma- 
trix into the Cartesian system T\, r 2 , r 3 , i.e., 

/(£>! + 2D 2 )/3 (£»i-£> a )/3 (D 1 - £> 2 )/3 \ 
T= (D!-D 2 )/3 (Di+2D 2 )/3 (Di - £> 2 )/3 . 
\(£)i-D 2 )/3 (Di-£> 2 )/3 (D 1 +2D 2 )/3J 

(20) 

This procedure requires the implementation of a full 
width matrix in the numerical evaluation of the frozen 
Gaussian thermal operator but avoids a full mass matrix 
and a change in the structure of the propagator. The 
results for the argon trimer presented in Sec. IIIIBI will 
show that despite its simplicity this choice leads to results 
which are competitive with the SP-TG propagator even 
though this variant of the frozen Gaussian propagator is 
much cheaper to evaluate numerically. 

To distinguish the frozen Gaussian propagator with the 
matrix structure of Eq. (120[) from its diagonal variant 
we will refer in the following to the two approximations 
as 2P-FG (two-parameter frozen Gaussian) and 1PD-FG 
(one-parameter diagonal frozen Gaussian), respectively. 



III. THE ARGON TRIMER 

A. Atomic parameters and numerical procedure 

To be able to compare our results with previous inves- 
tigations of the argon trimer^i^ we express the pairwise 
interaction by means of a Morse potential 

V(r t3 ) = D (exp [-2a(r l3 - R e )} - 2 exp [-a(r l3 - R e )]) 

(21) 

where is the distance between particles i and j. The 
Morse parameters are listed in Ref. [l6| and have been 
chosen such that the Morse potential reflects a previous 
fit to experimental results 2 ^. They are D = 99.00 cm -1 , 
a = 1.717 A, and R c = 3.757 A. Using the three sets 
of Gaussian parameters listed in Table U we achieved a 
very accurate Gaussian fit to this Morse potential with a 
standard deviation smaller than 0.4 cm -1 in the relevant 
region between r = 3.2 A and r — 6 A around the mini- 
mum. As will become clear below this deviation is much 
smaller than effects due to the Gaussian approximation 
used to calculate the quantum mean energy, i.e., differ- 
ences with the previous studies of the syste:m&i£ do not 
originate from the Gaussian fit of the potential, which is 
only introduced to accelerate numerical computations. 



In numerical evaluations of the partition function it 
is essential to restrict the position space integrations [cf. 
q integrations in Eqs. (fTUj) and (fl"5)) ] to a reasonable 
region of the configuration space containing all relevant 
information about the thermodynamics of the cluster. 
Usually, this is done by introducing an additional steep 
potential located at a certain distance R c from the center 
of mass -Rem, e.g., 

N , „ v 20 

^"Erx 1 (22) 

i=l ^ c ' 

(cf. Ref. 13) or, as implemented in our numerics, by re- 
stricting the sampling points q(r — 0) in the integrations 
to values \q — R cm \ < Rc- If the radius of the confining 
sphere is chosen correctly, the restriction or its explicit 
form should have no influence on the results. The Monte 
Carlo sampling in q is done with a standard Metropolis 
algorithm, where we followed the procedure suggested by 
Frantsuzov et al. explained in detail in Ref. 0. 

The optimum width parameter connected with the 
atom-atom interaction in the two parameter ansatz (f2"0")) 
of the frozen Gaussian propagator was found to be D 2 — 
25 A by using several choices, monitoring the results, 
and comparing them to the thawed Gaussian methods. 
For this choice of the parameter D 2 the low-temperature 
mean energy reaches the smallest value, i.e., a minimum 
when plotted vs. D 2 , thus, representing the best approxi- 
mation to the ground level. The observation of the small- 
est energy value in the limit T — > can be used as an 
additional criterion independently of the availability of 
a second method such as the thawed Gaussian propa- 
gator. We note that one can also monitor the relative 
amplitude of higher order corrections to the Gaussian 
approximatio n 15 ' 21 ' 24 as a function of the width param- 
eters. 

A value of D\ = 0.1 A 2 is already small enough to 
lead to the best possible description of the free center 
of mass motion. Using even smaller values for D\ did 
not change the results, so we decided to use this value to 
have a well conditioned matrix of which the eigenvalues 
do not differ too much in their order of magnitude. The 
parameter for a diagonal matrix in the 1PD-FG propa- 
gator representing the best middle ground between Di 
and D 2 is V = 20 A~ . It was selected with the method 
described for D 2 above. In a previous study of the frozen 
Gaussian imaginary time propagator it was found that 
the optimum choice for the width parameter is almost in- 
dependent of the temperature 1 -, this was also confirmed 
in our study of the argon trimer. We found that the ra- 
dius of the confining sphere does not have a significant 
influence on the optimum choice for the width matrices. 
Indeed, by checking several values as described above, 
it turned out that working with the same values for all 
computations was the best choice. 
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FIG. 1. Mean energy (a) and specific heat (b) of the argon 
trimer calculated for a confining radius of R c = 10 A. Results 
are provided for the thawed Gaussian approximations FC-TG 
and SP-TG, the frozen Gaussian approximations IPD-FG, 
2P-FG, and the classical theory. 



B. Comparison of the Gaussian imaginary time 
propagators 

First we investigated the different methods introduced 
in Sec. [IXJ for the evaluation of the quantum partition 
function for the argon trimer. To compare our results 
with the previous path-integral Monte-Carlo calculation 
by Perez de Tudela et al£ we selected one of their param- 
eter sets and calculated the mean energy and the specific 
heat for the argon trimer enclosed by a confining sphere 
with a radius of 10 A, which is the weakest confinement 
applied in their study. Our results obtained with the 
four Gaussian propagators described above are also com- 
pared with the corresponding derivatives of the classical 
partition function 



f kT \ 3/2N 
{ 2irh 2 ) 



3 -/JV(,) dq 3N 



(23) 



They are presented in Fig.[TJ To allow for the comparison 
with the path-integral Monte Carlo computations of Ref. 
6 in Fig. [5] we subtracted for this figure the exact kinetic 
energy of the free center of mass E cnl = 3/2kT from 
our values, which, in the calculation, always include the 
energy of the whole cluster including the center of mass 
translation. 

The four semiclassical methods are in reasonably good 
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FIG. 2. Comparison of the mean energy obtained with the 
Gaussian methods FC-TG, SP-TG, IPD-FG, and 2P-FG with 
the path-integral Monte-Carlo (PIMC) values taken from Ref. 
6. The kinetic energy of the free center of mass motion is 
subtracted from our results to allow for the comparison. 



agreement with each other and the classical results, how- 
ever, there are quantitative differences. One can ex- 
pect that the thawed Gaussian imaginary time propa- 
gator with a full matrix (FC-TG) provides the best ap- 
proximation to the exact quantum results, and indeed 
the mean energy obtained with that method shows the 
best correspondence with the path-integral Monte Carlo 
computations of Ref. (cf. Fig. [2]). In particular, in 
the low-temperature limit the FC-TG propagator gives 
the best approximation (E = — 246 cm -1 ) to the ground 
state energj£ of -252.44 cm" 1 . Even though the FC-TG 
method is the most versatile propagator applied here, it 
still assumes a Gaussian form of the atom's wave func- 
tions and shows a slight difference as compared to the nu- 
merically exact path-integral Monte Carlo ground state 
energy. For higher temperatures there appear larger de- 
viations as compared to the path-integral Monte Carlo 
results of Ref. [f| however, the latter ones are not very 
accurate due to the larger statistical error with increas- 
ing temperature. The semiclassical Gaussian approxima- 
tions are expected to describe the partition function and 
its derivatives at higher temperatures even better than at 
low temperatures since they converge to the correct an- 
swer in the classical limit. We can thus assume that our 
results are more accurate as the temperature increases. 
The differences between our results and those of Ref. H at 
higher temperatures T ^ 30 K become even clearer when 
considering the specific heat. Our Gaussian calculations 
indicate a direct transition from the bound cluster to the 
completely dissociated situation with three free atoms, 
whereas such a clear conclusion is not possible with the 
path-integral Monte Carlo calculations of Ref. [@, This 
transition is discussed in more detail in the next two sec- 
tions. 

It is also expected that the worst approximation in the 
low-temperature limit is obtained by the IPD-FG ansatz. 
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As was mentioned in Sec. Ill Bl the diagonal width matrix 
does not treat the free center of mass motion correctly. 
This manifests itself as a large deviation of the mean en- 
ergy from the correct value for T — > 0. The energy calcu- 
lated with the diagonal frozen Gaussian width matrix in- 
creases for temperatures below 5K, and this is definitely 
wrong. This deficiency is overcome with the 2P-FG ma- 
trix as presented in Eq. (I2TJ1) . It leads to a considerably 
lower mean energy in the low-temperature limit, which 
is closer to the exact ground level and is also closer to 
the FC-TG energy values, which can be considered to 
provide the best values of all methods used here. 

More interesting, however, is the comparison of the 
2P-FG and the SP-TG propagators. As reported in pre- 
vious investigation of larger clusters^ the single-particle 
thawed Gaussian approximation shows results for the 
mean energy and the specific heat which are qualitatively 
in agreement with the full matrix case. Quantitative dif- 
ferences have been reported and can also be found in our 
calculations for temperatures T J$ 45 K. The differences 
with respect to the full matrix results increase with de- 
creasing temperature. However, as can be seen in Fig. 
HJa) the description of the mean energy of the SP-TG 
and the 2P-FG methods is of similar quality at very low 
temperatures. This is a remarkable finding since the eval- 
uation of the single-particle thawed Gaussian propagator 
is more expensive than the frozen Gaussian due to the 
need to propagate the width matrix elements in time. 

By contrast, the specific heat curve of the 2P-FG prop- 
agator is closer to that of the 1PD-FG than to the two 
thawed Gaussian approximations which agree very well 
with each other. It seems that in the temperature re- 
gion around the transition from the bound to the disso- 
ciated cluster, a thawed Gaussian can provide a better 
description. However, the deviation is small and van- 
ishes for weaker external confinements. An example for 
such a weaker confinement is shown in Fig. [3J in which 
the comparison of Fig. Q] is repeated for a confinement of 
R c = 32 A. Here, down to a temperature of 18 K the two 
thawed Gaussian approximations and the 2P-FG propa- 
gator provide almost identical results. The lines lie on 
top of each other both for the mean energy and the spe- 
cific heat. Only the 1PD-FG values deviate a bit from 
the three other methods. 

In summary, we can state that the thawed Gaussian 
propagator with a full matrix provides the best approxi- 
mation which deviates in the low-temperature limit only 
slightly from the exact result and almost reaches the 
ground level for T — >• K. In all cases in which the high 
accuracy of the full matrix thawed Gaussian propagator 
is not required or in which the integration of the equa- 
tions of motion for the width parameters of a full matrix 
is too expensive, the frozen Gaussian ansatz with two pa- 
rameters (2P-FG) seems to be the best choice. It provides 
the same quality of results as the single-particle thawed 
Gaussian propagator but is much easier to evaluate since 
no integrations of parameters of the width matrix are 
required. 
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FIG. 3. Mean energy (a) and specific heat (b) of the argon 
trimer calculated with the Gaussian approximations FC-TG, 
SP-TG, 1PD-FG, 2P-FG, and classically for R c = 32 A. 



C. Dissociation and the influence of the confining sphere 

Before discussing the dissociation process in more de- 
tail we have to investigate the influence of the confining 
sphere on the mean energy and the specific heat. While 
often a very restrictive value of the confining radius R c 
is choseni^^^, our calculations demonstrate that its 
value significantly affects the thermodynamic properties 
of the clusters unless large values, drastically exceeding 
a few Angstroms, are used. Figure |4] shows the mean en- 
ergy of the argon trimer for several confinements R c in a 
purely classical calculation [Fig. 3(a)] and an evaluation 
of the quantum partition function with the 2P-FG prop- 
agator [Fig. H(b)], which, as was discussed in Sec. IIIIB1 
describes the quantum behavior correctly. The values 
R c = 2.6 A, 4 A, and 10 A have already been used in Ref. 
6. The qualitative behavior of these three curves differs 
strongly. In particular, the mean energy for R c = 10 A 
shows a significant increase of the slope for temperatures 
above 20 K. In Ref. [y this was regarded as an indication 
that this radius allows for a total fragmentation of the 
cluster. However, the mean energy for higher temper- 
atures reveals that this is not fulfilled completely. We 
included the line for E = (9/2)kT in the figure, which 
corresponds to three free particles. The mean energy 
for the restriction R c = 10 A does not reach that line 
even at T = 70 K. However, for the weaker restrictions 
R c = 15... 35 A one can observe that actually a total 
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FIG. 4. Comparison of the mean energy of the argon trimer 
for different confining radii R c (in A) in classical (a) and quan- 
tum calculations with the 2P-FG partition function (b). We 
also added the thin black lines E = 9/2kT representing three 
free particles and E = Vo + 6kT with the potential minimum 
Vo = —296 cm" 1 for the classical low-temperature behavior. 



fragmentation takes place and that for T gl 40 K the en- 
ergy is identical to that of three free particles. In Fig. 
SJa) we also included the line E — Vo + 6kT with the po- 
tential minimum Vq — — 296 cm -1 . This line corresponds 
to the classical expectation at low temperatures for the 
internal rotations and oscillations of the trimer plus the 
energy of the free center of mass. 

Even though the behavior at very low (T < 15 K) and 
high temperatures (T > 35 K) agrees well for all confine- 
ments R c > 15 A the transition itself obviously depends 
more critically on the value of R c . It is clear that a 
confining radius R c — 10 A is too restrictive and does 
not describe the cluster correctly. A significantly larger 
confining radius R c > 30 A is required. For the largest 
confining radii used in Fig. 0] we observe convergence, 
i.e., a further expansion of the confining sphere does not 
change the results significantly. 

We note that it becomes increasingly difficult to con- 
verge the Monte-Carlo integrations for the classical and 
the 2P-FG partition function for increasing R c . Sim- 
ilarly, the large error bars in the path-integral Monte 



Carlo calculation of Ref. HI indicate that in their com- 
putations a radius of R c — 10 A was already challenging. 
Nevertheless, the results presented in Fig. 2] demonstrate 
that the added effort of increasing R c beyond 30 A is 
essential. The necessity for a thorough investigation of 
the correct boundary conditions is already known from 
classical investigations of atomic clusters. Etters and 
Kaelberer— demonstrated the negative influence of too 
restrictive boxes on the classical average energy. 



D. Dissociation from classical and quantum perspectives 

The cluster at a confinement of R c = 32 A can be re- 
garded as converged with respect to R c . The artificial 
confinement does not have a further noticeable influence 
on the thermodynamic properties. This allows us to dis- 
cuss the features observed in the mean energy and the 
specific heat in more detail. Additionally, the various 
Gaussian propagators used to obtain the quantum prop- 
erties agree with each other to a high precision, so that 
we may consider them as converged in the sense that the 
choice of Gaussian method has no further influence. 

As can be seen in Fig.[3]the qualitative behavior in the 
quantum and classical cases is almost the same. The clus- 
ter is bound at low temperatures, shows a relatively sharp 
transition, and is completely dissociated for tempera- 
tures T > 33 K. For very low temperatures the classical 
mean energy exhibits the expected behavior E oc 6kT [cf. 
also Fig. Ufa)] for the system (free center of mass, rota- 
tions and oscillations of the internal degrees of freedom) , 
whereas the FC-TG mean energy (best approximation, 
see Sec. IIIIBI) converges to a value close to that of the 
ground level. Further differences between the classical 
and quantum mechanical results are found in the transi- 
tion region. It is shifted to slightly lower temperatures 
in the quantum calculations as compared to the classi- 
cal. All quantum calculations show a maximum of the 
specific heat at 20 K, while the classical maximum is at 
21.5 K. One reason for this shift is the presence of the 
zero point energy in the quantum system. As is already 
obvious in Fig. @] the classical and quantum results agree 
very well in the high-temperature limit, since both trend 
to the case of three free particles. Thus, we conclude 
that the transition observed in the cluster of three argon 
atoms is a classical phenomenon. The only difference be- 
tween classical and quantum mechanics is in the temper- 
ature at which the transition occurs. We note that Etters 
and Kaelberer— reported a "liquid-gas transition" iden- 
tified by the absence of bounded atom configurations in 
a classical investigation of the system with "free-surface 
boundary conditions", i.e., without a confining sphere, at 
T = 20 K, which is in good agreement with our results. 

Our finding of a classical-like complete dissociation of 
the trimer in one step differs from the conclusions of 
Perez de Tudela et al£. While the mean energy in the cal- 
culations of Ref. for R c = 10 A shows a larger and larger 
slope for increasing temperatures up to T — 40 K we ob- 
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serve already a decrease of the slope for temperatures 
T > 30 K. The difference relative to the path-integral 
Monte Carlo calculations becomes even more pronounced 
in the specific heat. Perez de Tudela et alA report that 
they find an "apparent" maximum which evolves with 
the radius R c of the confinement and appears slightly 
below 40 K for R c = 10 A. The absence of an unambigu- 
ous maximum was seen as an indication for structural 
changes of the cluster instead of a proper "phase transi- 
tion" . Although the dissociation of the cluster is not fully 
achieved for such a strong confinement, it is clear from 
the Gaussian methods used in this article that already for 
R c = I A a pronounced peak in the specific heat indicat- 
ing the dissociation of the system at T « 30 K is present 
(cf. Fig. Q]). Describing the cluster with a weaker con- 
finement correctly reveals the unambiguous dissociation 
without intermediate structural modification as discussed 
above. 

The very low temperature found for the dissociation 
of the cluster may also be important for rare gas clus- 
ters in general. As was already discussed^, even for a 
confining sphere with R c = 10 A, the transition temper- 
ature of T rj 35 K is lower than temperatures discussed 
for structural transformations or a "melting" of clusters. 
Features indicating such changes have, e.g., been found 
beyond 40 K for argon 25 . If one considers that the dis- 
sociation temperature is actually even lower (T w 20 K, 
cf. Fig. [3]), one necessarily concludes that it is wrong to 
ignore the influence of the confining sphere on such prop- 
erties. In larger neon clusters (Nei3 and Ne3s) features in 
the mean energy or the specific heat which were related to 
structural changes, have been reported between 6 K and 
8 Kli&I. These temperatures are lower than the dissocia- 
tion found here, however, one may expect that at least a 
partial dissociation can set in much earlier in larger clus- 
ters since they contain higher energies. In the numerical 
simulations^*^ the confining radii are chosen such that 
no atom can leave the cluster during the time evolution. 
Our results indicate that such a constraint might be too 
restrictive and lead to incorrect conclusions. A partial 
or full dissociation can influence structural transforma- 
tions and may even set in before structural changes of an 
artificially confined cluster can occur. 



IV. CONCLUSIONS AND OUTLOOK 

In the present article we investigated the argon trimer 
by means of semiclassical Gaussian approximations to 
the Boltzmann operator. We introduced a new matrix 
structure for a frozen Gaussian variant of the imaginary 
time propagator which is capable of correctly dealing 
with the free center of mass motion of a cluster of atoms 
in Cartesian coordinates. With this matrix structure we 
were able to show that the frozen Gaussian propagator 
is, in spite of its simplicity, competitive with numeri- 
cally more expensive thawed Gaussian variants. In par- 
ticular, the frozen Gaussian method provides the same 



quality thermodynamic results as the so-called single- 
particle thawed Gaussian propagator, which in addition 
to the time-dependent variables of the frozen Gaussian 
requires the time evolution for the elements of a block- 
diagonal width matrix. This is especially true in the low- 
temperature limit, where quantum effects become impor- 
tant and the form of the semiclassical approximation is 
supposed to have the largest influence. Our results sug- 
gest that the frozen Gaussian ansatz with two parameters 
(2P-FG) introduced in this article is the method of choice 
in all cases in which the higher accuracy of the full matrix 
thawed Gaussian propagator is not required or in which 
the integration of the equations of motion for the width 
parameters of a full matrix is too expensive. 



The evaluation of the mean energy and the specific 
heat for the cluster of three argon atoms demonstrated 
that a previous investigation of the system^ used too re- 
strictive confinements to describe the dissociation behav- 
ior of the system correctly. Above T = 15 K the cluster 
directly dissociates into three free atoms, as is evident 
from the Gaussian calculations presented in this article. 
This dissociation is almost purely classical, the influence 
of quantum mechanics is only found in a convergence to 
the ground state energy instead of the classical potential 
minimum for T — > and in a slight shift of the three- 
body dissociation transition to higher classical tempera- 
tures (AT Ri 1.5 K) which we attribute to the zero point 
energy of quantum mechanics. The clear and pronounced 
transition found in this article supports the conclusion of 
Ref. |6| that the dissociation of the atoms from the cluster 
is important when reconfigurations of the internal struc- 
ture are considered. Our results strongly indicate that 
the confinement to very small spheres usually applied in 
the calculation of the partition function and values de- 
duced from ikL&iil^ might be too restrictive to fully un- 
derstand the low-temperature behavior of the clusters. 
The dissociation can set in before structural changes or a 
melting can be observed. To make a clear statement on 
this question it is necessary to advance the investigations 
done here to clusters with higher numbers of atoms. In 
particular, the cases of Arg2&, Ari^^, Nei^ or Ne 38 i 
examined recently are of special interest. 



On the technical side, it is known that both the frozen 
and thawed Gaussian propagators used here are, in the 
framework of a generalized time-dependent perturbation 
theoryiaLSi, the lowest order approximations in a se- 
ries converging to the exact quantum propagato r 15 i 21 i 24 . 
Higher orders can help to understand the thermodynamic 
properties better and to verify the results obtained here 
with a higher accuracy. Furthermore, the corrections ob- 
tained by the evaluation of higher order terms provide 
objective access to the quality with which the Gaussian 
approximations reflect the quantum effects in the system 
studied. They are the topic of current studies. 
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